
library(hydroGOF)
library(clusterSim)
library(sqldf)

data=read.csv("C:/in/myd_merge.csv")
data=subset(data,data$station!=3)

data$aod=data.Normalization(data$aod,type="n5",normalization="column")
data$temp=data.Normalization(data$temp,type="n5",normalization="column")
data$avg_temp=data.Normalization(data$avg_temp,type="n5",normalization="column")
data$avg_rh=data.Normalization(data$avg_rh,type="n5",normalization="column")
data$avg_preci_24=data.Normalization(data$avg_preci_24,type="n5",normalization="column")

data1=data
model1=lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,data1)
#model1=lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,data1)
#model1=lm(data1$pm25~data1$aod+data1$temp+data1$avg_temp+data1$avg_rh+data1$avg_preci_24)

nSample=nrow(data1)
nPredict= length(model1$coefficients)-1

cutoff=4/(nSample-nPredict-1)
cook_thresold=qf(0.2,nPredict+1,nSample-nPredict-1)
cookValue=cooks.distance(model1)
data1["cookValue"] = cookValue

data2=subset(data1,data1$cookValue<=cook_thresold)

model2=lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,data2)
#model2=lm(data2$pm25~data2$aod+data2$temp+data2$avg_temp+data2$avg_rh+data2$avg_preci_24)

data2["predict"] = model2$fitted.values

pm25model = data2["predict"]
pm25real= data2["pm25"]
model_re=sum(abs(pm25real-pm25model)/pm25real)/nrow(data2)*100
model_rmse=rmse(pm25real,pm25model)
model_r=cor(pm25real,pm25model)

print("All Daily:")
print(paste("Number of sample: ",nrow(data2)))
print(paste("R2:",round(model_r*model_r,3)))
print(paste("RMSE:",round(model_rmse,3)))
print(paste("RE:",round(model_re,3)))

write.csv(data2,"C:/out/mod2.csv")

#---------------------------------------------------------------------

data=read.csv("C:/in/mod_daily_merge.csv")
data=subset(data,data$station!=3)

data$aod=data.Normalization(data$aod,type="n5",normalization="column")
data$temp=data.Normalization(data$temp,type="n5",normalization="column")
data$avg_temp=data.Normalization(data$avg_temp,type="n5",normalization="column")
data$avg_rh=data.Normalization(data$avg_rh,type="n5",normalization="column")
data$avg_preci_24=data.Normalization(data$avg_preci_24,type="n5",normalization="column")

data1=data
model1=lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,data1)
#model1=lm(data1$pm25~data1$aod+data1$temp+data1$avg_temp+data1$avg_rh+data1$avg_preci_24)

nSample=nrow(data1)
nPredict= length(model1$coefficients)-1

cutoff=4/(nSample-nPredict-1)
cook_thresold=qf(0.2,nPredict+1,nSample-nPredict-1)
cookValue=cooks.distance(model1)
data1["cookValue"] = cookValue

data3=subset(data1,data1$cookValue<=cutoff)

model3=lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,data3)
#model2=lm(data2$pm25~data2$aod+data2$temp+data2$avg_temp+data2$avg_rh+data2$avg_preci_24)

data3["predict"] = model3$fitted.values

#write.csv(data3,"C:/out/mod_daily.csv")


#---------------------------------------------------------------------

data=read.csv("C:/in/all_daily_merge.csv")
data=subset(data,data$station!=3)

data$aod=data.Normalization(data$aod,type="n5",normalization="column")
data$temp=data.Normalization(data$temp,type="n5",normalization="column")
data$avg_temp=data.Normalization(data$avg_temp,type="n5",normalization="column")
data$avg_rh=data.Normalization(data$avg_rh,type="n5",normalization="column")
data$avg_preci_24=data.Normalization(data$avg_preci_24,type="n5",normalization="column")

data1=data
model1=lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,data1)
#model1=lm(data1$pm25~data1$aod+data1$temp+data1$avg_temp+data1$avg_rh+data1$avg_preci_24)

nSample=nrow(data1)
nPredict= length(model1$coefficients)-1

cutoff=4/(nSample-nPredict-1)
cook_thresold=qf(0.2,nPredict+1,nSample-nPredict-1)
cookValue=cooks.distance(model1)
data1["cookValue"] = cookValue

data4=subset(data1,data1$cookValue<=cutoff)

model4=lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,data4)
#model2=lm(data2$pm25~data2$aod+data2$temp+data2$avg_temp+data2$avg_rh+data2$avg_preci_24)

data4["predict"] = model4$fitted.values
#write.csv(data3,"C:/out/all_daily.csv")

pm25model=model4$fitted.values
model_re=sum(abs(data4$pm25-pm25model)/data4$pm25)/nrow(data4)*100
model_rmse=rmse(data4$pm25,pm25model)
model_r=cor(data4$pm25,pm25model)

print("All Daily:")
print(paste("Number of sample: ",nrow(data4)))
print(paste("R2:",round(model_r*model_r,3)))
print(paste("RMSE:",round(model_rmse,3)))
print(paste("RE:",round(model_re,3)))
#------------------------------------------------------

modmyd_data=rbind(data2,data3)
#write.csv(modmyd_data,"C:/out/modmyd_daily.csv")
data_avg = sqldf("select station,aqstime,avg(predict) as avg_predict from modmyd_data group by station,aqstime;")
#write.csv(data_avg,"C:/out/modmyd_daily_group.csv")

#compare_data = sqldf("select data_avg.station,data_avg.aqstime,data_avg.avg_predict,data4.pm25 from data_avg inner join data4 on data_avg.station = data4.station and data_avg.aqstime = data4.aqstime;")
compare_data = sqldf("select data_avg.station,data_avg.aqstime,data_avg.avg_predict,data4.pm25,data4.predict from data_avg inner join data4 on data_avg.station = data4.station and data_avg.aqstime = data4.aqstime;")

model_r=cor(compare_data$avg_predict,compare_data$pm25)
model_re=sum(abs(compare_data$pm25-compare_data$avg_predict)/compare_data$pm25)/nrow(compare_data)*100
model_rmse=rmse(compare_data$pm25,compare_data$avg_predict)

print("MOD_MYD_daily:")
print(paste("Number of sample: ",nrow(compare_data)))
print(paste("R2:",round(model_r*model_r,3)))
print(paste("RMSE:",round(model_rmse,3)))
print(paste("RE:",round(model_re,3)))

model_r2=cor(compare_data$predict,compare_data$pm25)
model_re2=sum(abs(compare_data$pm25-compare_data$predict)/compare_data$pm25)/nrow(compare_data)*100
model_rmse2=rmse(compare_data$pm25,compare_data$predict)

print("All daily:")
print(paste("Number of sample: ",nrow(compare_data)))
print(paste("R2:",round(model_r2*model_r2,3)))
print(paste("RMSE:",round(model_rmse2,3)))
print(paste("RE:",round(model_re2,3)))





